function p_ri = ri(iter,beta,sigma,y,X,stations,tail,prob_treat)

    station1 = stations(:,1);
    station2 = stations(:,2);
    station3 = stations(:,3);
    
    sim_data = csvread('stations_randsim.csv',1,0);

    stationsim = sim_data(:,1);
    treatsim = sim_data(:,2:10001);

[N,k] = size(X);

treat1 = zeros(N, iter);
treat2 = zeros(N, iter);
treat3 = zeros(N, iter);

treatany = zeros(N,iter);
wgt=zeros(N,iter);
b = zeros(k+1,iter);
se=zeros(k+1,iter);
test = zeros(1,iter);

%Getting ac-wise treatment status
for j = 1:iter
    for i = 1:N
        for s = 1:60 
            if stationsim(s)== station1(i)
                treat1(i,j)=treatsim(s,j);
            end;
            if stationsim(s)==station2(i)
                treat2(i,j)=treatsim(s,j);
            end;
            if stationsim(s)==station3(i)
                treat3(i,j)=treatsim(s,j);         
            end;
            if treat1(i,j)+treat2(i,j)+treat3(i,j)>0
                treatany(i,j)=1;
            end;
        end;
    end;
end;

%Weights
for j = 1:iter
    for i = 1:N
        wgt(i,j)=(treatany(i,j))/prob_treat(i,1) + (1-treatany(i,j))/(1-prob_treat(i,1));
    end;
end;

for a = 1:iter
    [b(:,a), se(:,a)] = ipwreg(y,treatany(:,a),X,wgt(:,a),stations,tail);
end; 

%p-value
if tail == 1
    for l = 1:iter
        if b(1,l)/se(1,l) < beta(1)/sigma(1)
                test(1,l)=1;
          
        end;
    end;
    p_ri = mean(test,2);
end;

if tail == 2
     for l = 1:iter

             if abs(b(1,l)/se(1,l)) > abs(beta(1)/sigma(1))
                test(1,l)=1;
             end;
        
     end;
     p_ri = mean(test,2);
end;